#ifndef __MGGS_ELEMENT__
#define __MGGS_ELEMENT__

#include "Point.h"

#include <valarray>
#include <fstream>

template <typename T, int DIM>
class Element{
protected:
  int n_dofs;
  std::valarray<Point<T,DIM> > local_dofs;
  std::valarray<Point<T,DIM> > global_dofs;
  std::valarray<Point<T,DIM> > quad_points;
  std::valarray<double> quad_weights;
  int algebraic_accuracy;
public:
  Element();
  void set_n_dofs(const int _n);
  int get_n_dofs() const;
  void set_local_dofs(const int _i, const Point<T,DIM>& _P);
  const std::valarray<Point<T,DIM> >& get_local_dofs() const;
  const Point<T,DIM>& get_local_dofs(int _i) const;
  void set_global_dofs(const int _i, const Point<T,DIM>& _P);
  const std::valarray<Point<T,DIM> >& get_global_dofs() const;
  const Point<T,DIM>& get_global_dofs(int _i) const;
  int get_accuracy() const;
  int get_n_quad_points() const;
  const std::valarray<Point<T,DIM> >& get_quad_points() const;
  const Point<T,DIM>& get_quad_points(int _i) const;
  const std::valarray<double>& get_quad_weights() const;
  double get_quad_weights(int _i) const;
};

template <typename T, int DIM>
Element<T,DIM>::Element(){
  n_dofs = 0;
  algebraic_accuracy = 1;
}


template <typename T, int DIM>
void Element<T,DIM>::set_n_dofs(const int _n){
  n_dofs = _n;
  local_dofs.resize(_n);
  global_dofs.resize(_n);
}

template <typename T, int DIM>
int Element<T,DIM>::get_n_dofs() const{
  return n_dofs;
}

template <typename T, int DIM>
void Element<T,DIM>::set_local_dofs(const int _i, const Point<T,DIM>& _P){
  if (_i < n_dofs)
    local_dofs[_i] = _P;
  else
    std::cerr << "Valarray is out of bounds" << std::endl;
}

template <typename T, int DIM>
const std::valarray<Point<T,DIM> >& Element<T,DIM>::get_local_dofs() const{
  return local_dofs;
}

template <typename T, int DIM>
const Point<T,DIM>& Element<T,DIM>::get_local_dofs(int _i) const{
  return local_dofs[_i];
}

template <typename T, int DIM>
void Element<T,DIM>::set_global_dofs(const int _i, const Point<T,DIM>& _P){
  if (_i < n_dofs)
    global_dofs[_i] = _P;
  else
    std::cerr << "Valarray is out of bounds" << std::endl;
}

template <typename T, int DIM>
const std::valarray<Point<T,DIM> >& Element<T,DIM>::get_global_dofs() const{
  return global_dofs;
}

template <typename T, int DIM>
const Point<T,DIM>& Element<T,DIM>::get_global_dofs(int _i) const{
  return global_dofs[_i];
}


template <typename T, int DIM>
int Element<T,DIM>::get_accuracy() const{
  return algebraic_accuracy;
}

template <typename T, int DIM>
int Element<T,DIM>::get_n_quad_points() const{
  int n1 = quad_points.size();
  int n2 = quad_weights.size();
  if (n1 != n2){
    std::cerr << "Number of quad points error." << std::endl;
    std::exit(-1);
  }
  return n1;
}

template <typename T, int DIM>
const std::valarray<Point<T,DIM> >& Element<T,DIM>::get_quad_points() const{
  return quad_points;
}

template <typename T, int DIM>
const Point<T,DIM>& Element<T,DIM>::get_quad_points(int _i) const{
  return quad_points[_i];
}

template <typename T, int DIM>
const std::valarray<double>& Element<T,DIM>::get_quad_weights() const{
  return quad_weights;
}

template <typename T, int DIM>
double Element<T,DIM>::get_quad_weights(int _i) const{
  return quad_weights[_i];
}


template <typename T>
class TD_P1Element : public Element<T,2> {
private:
  int n_dofs = 3;
public:
  TD_P1Element();
  T local_basis_function(int _i, const Point<T,2> &_P) const;
  Vector<T,2>  local_basis_gradient(int _i) const;
  Point<T,2> local_to_global(const Point<T,2> &_P) const;
  Matrix<T,2,2> global_to_local_Jacobi() const;
  T global_to_local_det() const;
  Matrix<T,2,2> local_to_global_Jacobi() const;
  void read_quad_info(const char* filename, int _acc);
};

template <typename T>
TD_P1Element<T>::TD_P1Element():Element<T,2>(){
  this->Element<T,2>::set_n_dofs(n_dofs);
  Point<T,2> P0({0.0,0.0},-1,0);
  this->Element<T,2>::set_local_dofs(0,P0);
  Point<T,2> P1({1.0,0.0},-1,0);
  this->Element<T,2>::set_local_dofs(1,P1);
  Point<T,2> P2({0.0,1.0},-1,0);
  this->Element<T,2>::set_local_dofs(2,P2);
}

template <typename T>
T TD_P1Element<T>::local_basis_function(int _i, const Point<T,2> &_P) const{
  T xi = _P.get_coord(0);
  T eta = _P.get_coord(1);
  switch(_i){
    case 0:
	return (1 - xi - eta);
    case 1:
	return xi;
    case 2:
	return eta;
  }
}

template <typename T>
Vector<T,2>  TD_P1Element<T>::local_basis_gradient(int _i) const{
  switch(_i){
    case 0:
      return Vector<T,2>(-1.0,-1.0);
    case 1:
      return Vector<T,2>(1.0, 0.0);
    case 2:
      return Vector<T,2>(0.0, 1.0);
  }
}

template <typename T>
Point<T,2> TD_P1Element<T>::local_to_global(const Point<T,2> &_P) const{
  T x = 0;
  T y = 0;
  const std::valarray<Point<T,2> > &gdofs = Element<T,2>::get_global_dofs();
  for (int i = 0; i < n_dofs; i++){
    x += gdofs[i].get_coord(0) * local_basis_function(i, _P);
    y += gdofs[i].get_coord(1) * local_basis_function(i, _P);
  }
  Point<T,2> global_P({x,y},-1,0);
  return global_P;
}

template <typename T>
Matrix<T,2,2> TD_P1Element<T>::global_to_local_Jacobi() const{
  const std::valarray<Point<T,2> > &gdofs = Element<T,2>::get_global_dofs();
  Matrix<T,2,2> result;
  result(0, 0) = gdofs[1].get_coord(0) - gdofs[0].get_coord(0);
  result(0, 1) = gdofs[1].get_coord(1) - gdofs[0].get_coord(1);
  result(1, 0) = gdofs[2].get_coord(0) - gdofs[0].get_coord(0);
  result(1, 1) = gdofs[2].get_coord(1) - gdofs[0].get_coord(1);
  return result;
}

template <typename T>
T TD_P1Element<T>::global_to_local_det() const{
  Matrix<T,2,2> J = this->global_to_local_Jacobi();
  return J(0, 0) * J (1, 1) - J(0, 1) * J(1, 0);
}

template <typename T>
Matrix<T,2,2> TD_P1Element<T>::local_to_global_Jacobi() const{
  Matrix<T,2,2> result(2, 2);
  Matrix<T,2,2> J = global_to_local_Jacobi();
  T det = J(0, 0) * J (1, 1) - J(0, 1) * J(1, 0); 
  result(0, 0) = J(1, 1);
  result(0, 1) = -J(0, 1);
  result(1, 0) = -J(1, 0);
  result(1, 1) = J(0, 0);
  result /= det;
  return result;
}

template <typename T>
void TD_P1Element<T>::read_quad_info(const char* filename, int _acc){
  std::ifstream is(filename);
  int acc_table,acc;
  char garbage[64];
  is >> acc_table;
  if (acc_table < _acc){
    std::cerr << "The required algebraic accuracy " << _acc << " is too high,"
	      << " and the actual algebraic accuracy is " << acc_table << std::endl;
    acc = acc_table;
  }
  else
    acc = _acc;
  int next_acc,n_points;
  is >> next_acc;
  while (next_acc < acc){
    is >> n_points;
    for (int i = 0; i < n_points; i++)
      is.getline(garbage, 64);
    is.getline(garbage, 64);
    is >> next_acc;
  }
  is >> n_points;
  Element<T,2>::quad_points.resize(n_points);
  Element<T,2>::quad_weights.resize(n_points);
  for (int i = 0; i < n_points; i++){
    T x, y, w;
    is >> x >> y >> w;
    Element<T,2>::quad_points[i].set_coord({x,y});
    Element<T,2>::quad_weights[i] = w;
  }
  Element<T,2>::algebraic_accuracy = acc;
  is.close();
}


#else
// DO NOTHING.
#endif
